A criterion for the choice of optimal softening length $\epsilon$ and time-step $dt$ for $N$-body simulations of a collisionless stellar system is analyzed. Plummer and Hernquist spheres are used as models to follow how changes in various parameters of an initially equilibrium stable model depend on $\epsilon$ and $dt$. These dependences are used to derive a criterion for choosing $\epsilon$ and $dt$. The resulting criterion is compared to Merritt's criterion for choosing the softening length, which is based on minimizing the mean irregular force acting on a particle with unit mass. Our criterion for choosing $\epsilon$ and $dt$ indicate that $\epsilon$ must be a factor of 1.5-2 smaller than the mean distance between particles in the densest regions to be resolved. The time-step must always be adjusted to the chosen $\epsilon$ (the particle must, on average, travel a distance smaller than $0.5\epsilon$ during one time-step). An algorithm for solving N-body problems with adaptive variations of the softening length is discussed in connection with the task of choosing $\epsilon$, but is found not to be promising.